A three-dimensional wavelet based multifractal method : about the need of revisiting 
the multifractal description of turbulence dissipation data 
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We generalize the wavelet transform modulus maxima (WTMM) method to multifractal analysis 
of 3D random fields. This method is calibrated on synthetic 3D monofractal fractional Brownian 
fields and on 3D multifractal singular cascade measures as well as their random function counterpart 
obtained by fractional integration. Then we apply the 3D WTMM method to the dissipation field 
issued from 3D isotropic turbulence simulations. We comment on the need to revisiting previous 
box-counting analysis which have failed to estimate correctly the corresponding multifractal spectra 
because of their intrinsic inability to master non-conservative singular cascade measures. 
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The multifractal formalism was introduced in the mid- 
eighties to provide a statistical description of the fluctu- 
ations of regularity of singular measures found in chaotic 
dynamical systems 1] or in modelling the energy cascad- 
ing process in turbulent flows Q, IH HI- Box-counting 
(BC) algorithms were successfully adapted to resolve 
multifractal scaling for isotropic self-similar fractals |^|^. 
As to self-affine fractals, Parisi and Frisch developed, 
in the context of turbulence velocity data analysis, an al- 
ternative multifractal description based on the so-called 
structure functions. Unfortunately, there are some draw- 
backs to these classical multifractal methods and as pro- 
posed in Ref. a natural way of performing a unified 
multifractal analysis of both singular measures and multi- 
affine functions, consists in using the continuous wavelet 
transform (WT). Applications of the so-called WTMM 
method to ID signals have already provided insight into 
a wide variety of problems, e.g. fully-developed turbu- 
lence, financial markets, meteorology, physiology, DNA 
sequences Q. Recently, the WTMM method has been 
generalized to 2D for multifractal analysis of rough sur- 
faces 0, and successfully applied to characterize the in- 
termittent nature of satellite images of the cloud struc- 
ture 'lO'l and to assist in the diagnosis in digitized mam- 
mograms . Our aim here is to go one step further and 
to generalize the WTMM method from 2D to 3D. 

The 3D WTMM method consists in smoothing the dis- 
crete 3D field data by convolving it with a filter and 
then in computing the gradient on the smoothed sig- 
nal as for multiscale edge detectors j^^. Define three 
wavelets %lji{x,y,z) = d4>{x,y, z)/dxi with Xi = x,y or z 
for i = 1, 2 or 3 respectively and 4'{x, y, z) is a 3D smooth- 
ing function well localized around x = y — z — Q. For 
any function f{x,y, z) £ L^(M''), the WT at the point b 
and scale a can be expressed in a vectorial form Il2l| : 

T^[/](b,a) = V{0b,a*/}, (1) 

where 0b,a(r) — a^^0(a^-^(r — b)). If tj) is just a Gaus- 
sian 0(r) = exp(— r^/2) (or one of its derivatives), then 
Eq. defines the 3D WT as the gradient field vector of 



/(r) smoothed by dilated versions (p{r/a) of this filter. At 
a given scale a, the WTMM are defined by the positions 
b where the WT modulus 7W^[/](b,a) = |T^[/](b,a)| 
is locally maximum along the direction of the WT vec- 
tor (Eq. (P)). These WTMM lie on connected surfaces 
called maxima surfaces (see Fig. 12}. In theory, at each 
scale a, one only needs to record the position of the local 
maxima of (WTMMM) along the maxima surfaces 
together with the value of and the WT vector 

direction. They indicate locally the direction where the 
signal has the sharpest variation. These WTMMM are 
disposed along connected curves across scales called max- 
ima lines 9, 12] living in a 4D space (x, y, z, a). We will 
define the WT skeleton as the set of maxima lines that 
converge to the (x, y, z ) -hyperplane in the limit a — > 0^ 
(FiglHd)). One can prove [sUlJI that, provided the first 
n.0 moments off/? be zero, then A4,p[f][Cro{a)^ ^ a''^""") 
along the maxima line £ro pointing to the point Tq in 
the limit a 0"*", where h{ro) (< JJ^iis the local Holder 
exponent of /. As in ID and 2D the 3D WTMM 

method consists in defining the partition functions : 

Z{q,a)= J2 (A^^[/](r,«))'^a-(^), (2) 

£6£(a) 

where q £ W and C{a) is the set of maxima lines of the 
WT skeleton. Then from Legendre transforming T(q): 
D{h) — muiq{qh — T{q)^, one gets the D{h) singularity 
spectrum defined as the Hausdorff dimension of the set of 
points r where h{r) is h. As an alternative strategy one 
can compute the mean quantities h(q, a) and D{q, a) : 

h{q,a)= ln|M^[/](r,a)| W^^[/](g,Aa) , (3) 

CeC(a) 

D{q,a)^ W^[f]{q,C,a)\n{W^[f]{q,C,a)) , {A) 

where W.4,[f]{q, C, a) = {M^|,[f]{r, a)y / Z{q, a) is a 
Boltzmann weight computed from the WT skeleton. 
From the scaling behavior of these quantities, one 
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can extract h{q) = lima^Q+ h{q, a)/ \n a and D{q) = 
lmia^o+ D{q,a)/lna and therefore the D(h) spectrum. 

Fractional Brownian motions (fBm) are homogenous 
random self-affine functions that have been used to cal- 
ibrate both the ID and 2D WTMM methods 3D 
fBm Bh[y) are Gaussian stochastic processes with sta- 
tionary increments and well-known statistical properties 
: r(g) = qH — 3, 0<iJ<l. In Fig. 1 are reported the 
results of the 3D WTMM method when applied to 16 x 
(256)'^ realizations of Bh=i/2{^)- As shown in Fig. 1(a), 
Z{q,a) (Eq. (O) display nice scahng behavior over 3 oc- 
taves when plotted t^s a in a logarithmic representation. 
A linear regression fit of the data yields the linear r(g) 
spectrum shown in Fig. Efc), in good agreement with the 
theoretical spectrum. This signature of monofractality is 
confirmed in Fig. ^b) where, when plotting h{q, a) vs 
log2(a), H — 1/2 is shown to provide an excellent fit of 
the slope of the data {aw ^ a ^ 4:aw) for q e] — 2,4[. 
When using h{q,a) (Eq. (0)) and D{q,a) (Eq. (g))) to 
estimate h{q) and D{q) in Fig. G^d), one gets, up to the 
numerical uncertainty, a D(h) spectrum that reduces to 
a single point D(h = H = 1/2) = 3. Similar quantitative 
estimates have been obtained for H = 1/3 and 2/3, thus 
confirming that 3D fBm's are nowhere differentiable with 
a unique Holder exponent h = H. 

Generating multifractal measures using multiplicative 
cascades is well documented ^;2;i, ;.3{, :J] . The "binomial (or 
p-) model", originally designed to account for the sta- 
tistical scaling properties of the dissipation field in fully 
developed turbulence, has very simple multifractal prop- 
erties Q 0, 0] ■ The 3D version of the p-model consists 
in starting with a cube of size L, in which a measure 
/i = /ii is uniformly distributed. At the first step, the 
initial cube is broken into eight smaller cubes of size L/2 
and one selects at random the four sub-cubes which will 
receive a fraction M^^^ = pi/4, the four others receiving 
the fraction M^^^ = p2/4 of the measure (pi + P2 — 1)- 
Iterating this rule, one generates a random singular mea- 
sure Ai„(r;0 = t^LU^=lM^ , l/L = 2-« ^ 0. A 
straightforward computation yields : Tf^{q) — ~{q + 
2) — log2(pf -I- pI)- As a mean of introducing continu- 
ity, Fractionaly Integrated Singular Cascade (FISC) al- 
gorithm [l^ amounts to a low-pass power-law filtering in 
Fourier space: /„(r) = ^„(r) * |r|"(^~^') , < H* < 1. 
This leads to FISC random functions with the follow- 
ing multifractal spectrum : T/(g) = T^{q) + qH* = 
-2 - g(l - H*) - log^ipl+pl). In Fig. fflare reported 
the results of the 3D WTMM analysis of the FISC model 
with pi = 0.32 and H* = 0.638. As shown in Fig. da), 
Z{q, a) display good scaling for q g] — 2, 4[ for which sta- 
tistical convergence turns out to be achieved. The corre- 
sponding Tf{q) spectrum is displayed in Fig. ^c) along 
with the theoretical spectrum. The agreement is quite 
satisfactory; Tf(q) is nonlinear, the hallmark of multi- 
fractal scaling. This is confirmed in Fig. ^b) where the 
slope of h{q, a) vs log2(a) clearly depends on q. From the 
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FIG. 1: Multifractal analysis of 3D (■)! 3D p-model 

(•, p = 0.32) and 3D FISC (o, p = 0.32, H* = 0.638) data 
using the 3D WTMM method with a first (■) or third (•,o) 
order analyzing wavelet (</>(r) is a Gaussian or its Laplacian). 
(a) logj Z[q, a) vs logj a; (b) h[q, a) vs logj a; the dashed and 
solid lines correspond to linear regression fits over the three 
first octaves, (c) r(g) vs q; the dashed and solid lines are the 
theoretical predictions, (d) D{h) vs h as obtained from the 
scaling behavior of h{q,a) and D{q,a). These results corre- 
spond to annealed averaging over 16 x (256)^ images, a is 
expressed in aw units where aw ~ 7 (pixels) is the charac- 
teristic size of at the smallest resolved scale. In (c) and (d) , 
the symbols (a) represent the results of classical box-counting 
analysis of the 3D p-model. 



estimate of h{q) and D{q), one gets the single- humped 
Df{h) spectrum shown in Fig. Q^d). Note that some de- 
parture from the theoretical spectrum can be observed on 
the right-hand side of the D(/i)-curve {q < —2) indicat- 
ing that the estimate of the weakest singularities would 
require a larger statistical sample. In Figs, ^c) and^d) 
are also reported the results of a comparative analysis of 
the p-model using the 3D WTMM method and classical 
BC techniques. Note that the WTMM definition of the 
T^^(g) spectrum (Eq. @) shghtly differs from the BC 
definition T^*^(g) found in the literature I3| : 



T^''{q)^T::'"{q) + dq={q-m 
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which implies the following relationship ff^'^{a) = 
D'^'^ {h = a — d) between the corresponding singular- 
ity spectra {d — 3). For both methods, the numerical 
results for T^{q) and D^{h) are found in good agreement 
with the theoretical spectra, for q €] — 2,4[. In particu- 
lar, the cancellation exponent 14] is found t^'^((? = 1) = 
rj^^iq = l)-l-3 = 0, as the signature of the conservativity 
of the p-model cascading rule. One of the main problem 
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with the BC method is the fact that, by construction, 
the measure in a given box is the sum of the measures 
in smaller non-overlapping boxes, which implies that the 
cancellation exponent t^*^((7 = 1) = 0. This means that 
BC algorithms are not adapted to study non-conservative 
singular cascades, signed measures as well as multifractal 
functions for which the cancellation exponent has no rea- 
son to vanish. Altogether the results reported in Fig. 1 
bring the demonstration that our 3D WTMM methodol- 
ogy paves the way from multifractal analysis of singular 
measures to continuous multi-affine functions. In partic- 
ular the Df{h) curve for the FISC random function is 
found identical to the D^{h) curve up to a translation to 
the right by H* ; this is the consequence of the fractional 
integration which implies Df{h) = D^{h — H*). 

Since Kolmogorov's founding work (K41), fully devel- 
oped turbulence has been intensively studied theoreti- 
cally, numerically and experimentally 0, A central 
quantity in the K41 theory is the mean energy dissipation 
e = J2i which is supposed to be constant. 

Indeed, e is not spatially homogenous but undergoes lo- 
cal intermittent fluctuations [^^. There have been early 
experimental attemps to measure the multifractal spec- 
tra of e 13 ■ Surprisingly, the binomial model turns out to 
account reasonably well for the observed (q) and (a) 
spectra Q ■ Experimentally, single probe measurement of 
the longitudinal velocity requires the use of the ID sur- 
rogate dissipation approximation e' = 15h'{du/dx)^ that 
may introduce severe bias in the multifractal analysis. 
3D multifractal processing of dissipation data is at the 
moment feasible only for numerically simulated flows at 
moderate Reynolds number for which scaling just begins 
to manifest itself 0. Several numerical studies [16a, c] 
agree that e' is in general more intermittent than e which 
is found nearly log-normal in the inertial range [16b, c]. 

So far mainly BC techniques have been used to per- 
form multifractal a naly sis of numerical and experimental 
dissipation data [lElj- Here we apply the 3D WTMM 
method to isotropic turbulence DNS data obtained by 
Meneguzzi with the same numerical code as in Ref. M 
but at a (512)'^ resolution and a viscosity of 5.10"''' cor- 
responding to a Taylor Reynolds number Rx — 216 (one 
snapshot of the dissipation 3D spatial field). For the sake 
of comparison, we will also report the results of some av- 
eraging over 18 snapshots of (256)'^ DNS run by Leveque 
at R\ = 140. The main steps of our 3D WT computation 
are illustrated in Fig. |21 Focusing on a (64'^) sub-cube, 
we show the original e data (Fig. Eta)), the WTMM sur- 
faces along with the WTMMM points computed at two 
different scales (Figs. |2Ib) and|2Ic)) and some projection 
of the WT skeleton (Fig. E^d)). According to Eq. ©, by 
plotting Z{q, a)/{Z{q = 0, a)Y a in a logarithmic rep- 
resentation in Fig. El^a), one can then directly compare 
{Z[q = 0, a) - a^"^) our WTMM computations with BC 
ones. Actually, good scaling properties are observed for 
q s] — 2,4[. Linear regression fits of the data yield the 




FIG. 2: 3D WT analysis of the DNS dissipation field e. i/> is 
the first-order analyzing wavelet {4>{r) is the Gaussian), (a) 
Isosurface plot of e in a (64)^ sub-cube, (b) WTMM surfaces 
at scale a = 2(tw\ from the local maxima (WTMMM) of JVi^ 
along these surfaces originate a black segment whose length 
is proportional to and direction is along the WT vector, 
(c) same as in (b) for a = Aaw- (d) 3D projection of the WT 
skeleton obtained by linking the WTMMM across scales. 



non-linear T^{q) spectra shown in Fig. |3|[c) that signifi- 
cantly deviate from a straight line, the hallmark of niul- 
tifractality. But surprisingly, rj^^ (g) significantly differs 
from T^^{q) - 'dq. Actually, our 3D WTMM algorithms 
reveal that the cancellation exponent is significantly 
different from zero : rf''^ {q = 1) 3 = -0.19 ± 0.03 < 0, 
the signature of a signed measure. Indeed, as shown 
in Fig. Etc), T^^{q) data are rather nicely fitted by 
the theoretical spectrum T^{q) of the non-conservative 
p-model with pi — 0.36 and p2 — 0.78 {pi + P2 ~ 
1.14 > 1). Thus BC algorithms systematically provide 
a misleading conservative (g) spectrum diagnostic with 
P = Pi/ {Pi + P2) and I - p = P2I {p\ + P2)- As shown 
in Fig. Ofc), T^'~^(q) — 3q data are quite well reproduced 
by the theoretical conservative p-model spectrum with 
p = 0.32 = pi/{pi + P2) = 0.36/1.14, consistently with 
our 3D WTMM finding for r,^^(g). The difference be- 
tween the two spectra is nothing but a fractional integra- 
tion of exponent H* — log2(pi + P2) ^ 0.19. This result 
is confirmed in Fig. Old) where the singularity spectrum 
f^'~^ {a) is misleading shifted to the right by H* (= — the 
cancellation exponent) ~ 0.19, without any shape change 
as compared to the [a) = D^'^{h = a — d). Note 
that the observation that fj^'^{a — d + 1) is not tan- 
gent to the diagonal (on the contrary to f^'~''{a — d+1)) 
is clearly related to the fact that the cancellation expo- 
nent is different from zero [T^ . This observation seriously 
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questions the validity of most of the experimental and nu- 
merical BC estimates of the t^'~^ {q) and f^'~'{a) spectra 
reported so far in the literature. In Fig. |3td) is shown 
for comparison some average fF^ (a) spectrum obtained 
by Meneveau and Sreenivasan Q for ID surrogate dissi- 
pation experimental data {d = 1). We notice that these 
experimental data were claimed to be well fitted by the 
conservative p-model with p = 0.3 (1 — p = 0.7), i.e. a 
value not so far from the value p — pi/{pi + P2) = 0.32 
derived from our 3D WTMM method. 

Previous application of the ID WTMM method to ID 
surrogate dissipation data has already revealed the fact 
that the cancellation exponent might be different from 
zero [ill- The 3D WTMM resuhs reported here show 
that this surprising result is not some artefact resulting 
from the analysis of ID cuts (the dissipation could well 
be non conserved along these cuts), but that the mul- 
tifractal spatial structure of the 3D dissipation field is 
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FIG. 3: Multifractal analysis of Meneguzzi DNS dissipation 
data (d = 3) using the 3D WTMM method (A) and BC 
techniques (a). (a) \og2{Z{q,a) / {Z{q = 0, o))^) vs logj a; 

(b) h{q, a) — a{q, a) — 3 vs logj a; the solid and dashed lines 
correspond to linear regression fits over aw ^ o ^ 

(c) Ti:{q) = Tj*^(q) or r^^' (q) — 3q vs g; the solid (dashed) 
lines correspond to the p-model prediction with pi — 0.36, 
P2 = 0.78 (pi = 0.32, p2 = 0.68). (d) /,(a) -f 3 - d vs 
a-d+1, where /,(a) = /f'^(a) or /J^^(a) = DY^{h = a-d) 
; the solid and dashed lines have the same meaning as in 
(c); the thick solid line is the log-normal spectrum fe{a) — 
d - {a-d+1- Ci)V2C2 with Ci = 0.91 and C2 = 0.22. 
The symbols (□) correspond to some average /(a) spectrum 
of experimental (d = 1) surrogate dissipation data In (c) 
and (d), the 3D WTMM muhifractal spectra of Leveque 3D 
dissipation data (18 snapshots) (o) are shown for comparison. 
In (d) the dashed straight line is the diagonal. 

likely to be well described by a multiplicative cascade 



process that is definitely non conservative. As shown in 
Fig. Ofd) , this conclusion is confirmed by the results ob- 
tained when averaging over 18 snapshots of Leveque's 
DNS; the cancellation exponent is found even more neg- 
ative T^^{q — 1) -I- 3 = —0.26, as an indication (as re- 
gards to the smaller R\ value) that this exponent might 
decrease to zero in the limit of infinite Reynolds number. 
To conclude, let us point out that the f^'^{a) data seem 
to be even better fitted by a parabola, as predicted for 
non-conservative log-normal cascade processes. 

We are very grateful to M. Meneguzzi and E. Leveque 
for allowing us to have access to their DNS data and to 
the CNRS under GDR turbulence. 
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